#ifndef LAMINAR_FLOW_HH
#define LAMINAR_FLOW_HH

#include "laminarFlow.h"

#include <vtkVersion.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderWindow.h>
#include <vtkSmartPointer.h>
#include <vtkChartXY.h>
#include <vtkTable.h>
#include <vtkPlot.h>
#include <vtkFloatArray.h>
#include <vtkContextView.h>
#include <vtkContextScene.h>
#include <vtkPen.h>

#include <vtkAxis.h>
#include <vtkWindowToImageFilter.h>
#include <vtkPNGWriter.h>

namespace cfd {

// Constructor
template<typename T>
LaminarFlow<T>::LaminarFlow(int n, T mu_const, T L, T p_0, T p_L, T c) : 
  _n(n), _mu_const(mu_const), _L(L), _p_0(p_0), _p_L(p_L), _c(c), 
  _grad_p(-(p_L-p_0)/L),
  _A(n, vector<T>(3,0)), _w(n,0), _b(n,0),
  _mu(n+1, mu_const), _h(_n), _r(n+1,0),
  _it(0) 
{
  for (int i=1; i<(_n+1); i++) {
    _r[i] = 1.0*i/_n;
    _h[i-1] = _r[i]-_r[i-1];
  }
}

// Destructor
template<typename T>
LaminarFlow<T>::~LaminarFlow() {
}

// One iteration step
template<typename T>
void LaminarFlow<T>::oneStep() {
  int i; 
  
  _grad_p = -(_p_L-_p_0)/_L*cos(0.1*_it); // dp/dz
  _it++;

  // 2.1 contorno interior: i=0
  i = 0; // ignore first column
  _A[0][1] = (_mu[i+1]*_r[i+1]) / (_r[i+1]+0.5*_h[i+1]) + _r[i+1]*_r[i+1]*_c;
  _A[0][2] = -(_mu[i+1]*_r[i+1]) / (_r[i+1]+0.5*_h[i+1]);
  _b[0] = _grad_p*0.5*_r[i+1]*_r[i+1] + _r[i+1]*_r[i+1]*_c*_w[i];
  
  // 2.2 dentro: i = 1:(n-2)
  for (i=1; i<(_n-1); i++) {
    _A[i][0] = -(_mu[i]*_r[i]) / (0.5*(_h[i]+_h[i-1]));
    _A[i][2] = -(_mu[i+1]*_r[i+1]) / (0.5*(_h[i+1]+_h[i]));
    _A[i][1] = -_A[i][0] - _A[i][2] + (_r[i+1]*_r[i+1] - _r[i]*_r[i])*_c;
    _b[i] = _grad_p*0.5 * (_r[i+1]*_r[i+1]-_r[i]*_r[i]) + (_r[i+1]*_r[i+1] - _r[i]*_r[i])*_c*_w[i];
  }
  
  // 2.3 contorno externo: i = (n-1)
  i = _n-1;
  _A[i][0] = -(_mu[i]*_r[i]) / (0.5*(_h[i]+_h[i-1]));
  _A[i][1] = -_A[i][0] + 2* (_mu[i+1]*_r[i+1]) / _h[i] + (_r[i+1]*_r[i+1]-_r[i]*_r[i])*_c;
  _b[i] = _grad_p*0.5*(_r[i+1]*_r[i+1]-_r[i]*_r[i]) + (_r[i+1]*_r[i+1]-_r[i]*_r[i])*_c*_w[i];
  
  // 2.4 Solve Aw=b
  resolveLES();
  //showVector(&_w);
  
//  plot();
}

// destroys the diagonal of A, i.e. _A[i][1], and the right-hand side, i.e. _b[i]
template<typename T>
void LaminarFlow<T>::resolveLES() {
  // Thomas algorithm
  for (int i = 1; i < _n; i++) {
    double m = _A[i][0]/_A[i-1][1];
    _A[i][1] = _A[i][1] - m * _A[i-1][2];
    _b[i] = _b[i] - m*_b[i-1];
  }
 
  _w[_n-1] = _b[_n-1]/_A[_n-1][1];

  for (int i=_n-2; i>=0; --i)
    _w[i] = (_b[i]-_A[i][2] * _w[i+1]) / _A[i][1];  
}

// Plot one cycle
template<typename T>
void LaminarFlow<T>::plot(){
  // Create a table with some points in it
  vtkSmartPointer<vtkTable> table = vtkSmartPointer<vtkTable>::New();
  
  vtkSmartPointer<vtkFloatArray> arrZ = vtkSmartPointer<vtkFloatArray>::New();
  arrZ->SetName("z (distance to center)");
  table->AddColumn(arrZ);
 
  vtkSmartPointer<vtkFloatArray> arrW = vtkSmartPointer<vtkFloatArray>::New();
  arrW->SetName("w (velocity)");
  table->AddColumn(arrW);
 
  // Fill in the table with some example values
  table->SetNumberOfRows(_n);

  oneStep();
  table->SetValue(0, 0, 0.0);
  table->SetValue(0, 1, _w[0]);

  for (int i = 1; i<_n; ++i) {
    table->SetValue(i, 0, _r[i]+0.5*_h[i]);
    table->SetValue(i, 1, _w[i]);
  }

  // Set up the view
  vtkSmartPointer<vtkContextView> view = vtkSmartPointer<vtkContextView>::New();
  view->GetRenderer()->SetBackground(1.0, 1.0, 1.0);
  
  // Add multiple line plots, setting the colors etc
  vtkSmartPointer<vtkChartXY> chart = vtkSmartPointer<vtkChartXY>::New();
  view->GetScene()->AddItem(chart);

  chart->GetAxis(vtkAxis::LEFT)->SetTitle("w (velocity)");
  chart->GetAxis(vtkAxis::LEFT)->SetRange(-.25, .25);
  // Prevent rescaling the axis
  chart->GetAxis(vtkAxis::LEFT)->SetBehavior(vtkAxis::FIXED);
  chart->GetAxis(vtkAxis::BOTTOM)->SetTitle("z (distance to center)");
  chart->SetShowLegend(true);

  vtkPlot *line = chart->AddPlot(vtkChart::LINE);
  line->SetInput(table.GetPointer(), 0, 1);
  line->SetColor(0, 255, 0, 255);
  line->SetWidth(1.0);
  //line->GetPen()->SetLineType(2);//For dotted line, can be from 2 to 5 for different dot patterns
  
 
  view->GetRenderWindow()->SetMultiSamples(0);
 
  // Start interactor
  //view->GetInteractor()->Initialize();
 
  //view->GetInteractor()->Start();
  //view->Render();
  char buffer[256];

  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(view->GetRenderer());
  renderWindow->SetSize(320,240);
//  renderWindow->OffScreenRenderingOn(); 
//  printf("\n%8.4f\n",_w[_n-1]);
//  showVector(&_w);
  
  double tmp[2];
  tmp[0] = _w[0];
  int flipped = 0;
  
  for (int j=1; flipped<3; j++) {
//  for (int j=1; j<100; j++) {
    // create tables
    table->SetValue(0, 0, 0.0);
    table->SetValue(0, 1, _w[0]);

    for (int i = 1; i<_n; ++i) {
      table->SetValue(i, 0, _r[i]+0.5*_h[i]);
      table->SetValue(i, 1, _w[i]);
    }

    // Redraw plot
    chart->RemovePlot(vtkChart::LINE);
    line = chart->AddPlot(vtkChart::LINE);
    line->SetInput(table.GetPointer(), 0, 1);
    line->SetColor(0, 255, 0, 255);
    line->SetWidth(1.0);
    //view->Render();

    // write normal sized image
    // Render plot and save it to file
    if (flipped>0) {
      vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter = vtkSmartPointer<vtkWindowToImageFilter>::New();
      windowToImageFilter->SetInput(renderWindow);
      windowToImageFilter->Update();

      vtkSmartPointer<vtkPNGWriter> writer = vtkSmartPointer<vtkPNGWriter>::New();
      sprintf(buffer, "tmp/image_%03d.png", _it);

      writer->SetFileName(buffer);
      writer->SetInput(windowToImageFilter->GetOutput());
      writer->Write();
    }
    
    oneStep();
    tmp[j%2] = _w[0];
    if (tmp[0]*tmp[1]<0)
      flipped++;
    }

}

template<typename T>
void LaminarFlow<T>::showVector(vector<T> *vec){
  typename vector<T>::iterator i;
  for(i=vec->begin(); i!=vec->end(); i++) {
      printf("%8.5f\n", *i);   
  }
  printf("\n");
}

template<typename T>
void LaminarFlow<T>::printA() {
  typename vector< vector<T> >::iterator i;
  typename vector<T>::iterator j;
  for(i=_A.begin(); i!=_A.end(); i++) {
    for(j=(*i).begin(); j!=(*i).end(); j++) {
      printf("%8.5f, ", *j);
    }
    printf("\n");
  }      
}

} // namespace cfd


#endif // LAMINAR_FLOW_HH


